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We derive two sets of explicit algebraic constraint preserving boundary conditions for the lin- 
earized BSSN system. The approach can be generalized to inhomogeneous differential and evolution 
conditions, the examples of which are given. The proposed conditions are justified by an energy 
estimate on the original BSSN variables. 
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I. INTRODUCTION 

The widely used treatment of Einstein's equations 
in numerical relativity is to cast them to the form of 
a nonlinear hyperbolic system with constraints (e.g., 
01 HIS El El) and solve by employing sophisticated dis- 
cretization techniques. In the course of solution, the con- 
straint part is either monitored, or explicitly imposed. It 
was observed, that the solution of the evolution part with 
no constraints produces a violation which grows rapidly 
breaking computations in a short time 0-0 An attempt 
to control constraint violation, by projecting the solution, 
or by incorporating constraint quantities in the evolution 
equations, results in a longer life time of calculations as 
a rule of thumb (e.g., 0, HIS 01)- It was found in 
that exponentially growing constraint violating solutions 
converge to unstable solutions of the dynamic equations, 
which suggests that the constraint violation is closely re- 
lated to loss of stability in the system. 

An exact solution to evolution equations in the entire 
space has a property that it satisfies constraint equations 
automatically as long as it satisfies them initially. How- 
ever, in numerical simulations, because of the roundoff 
and truncation errors, one cannot hope for automatic 
constraint compliance. Instead, care must be taken to 
ensure that the inserted perturbations are small, and re- 
main small during the evolution. 

The behavior of the solution can be improved signif- 
icantly H [nj by introducing special sets of boundary 
data, the so-called constraint-preserving boundary con- 
ditions, or conditions that imply trivial evolution of con- 
straints. Several sets of such data were proposed for var- 
ious first order formulations of Einstein's equations (e.g., 
d [Hill ESQ) 1 These conditions are typically writ- 
ten as a system of partial differential equations restricted 
to the boundary, and in cases when the equations are 
time dependent and decouple from the bulk system, the 
equations may be integrated in time to produce regular 



Dirichlet data that is compatible with constraints [l2| . 

In this work, two sets of well-posed homogeneous alge- 
braic constraint-preserving boundary conditions for the 
linearized Baumgarte-Shapiro-Shibata-Nakamura for- 
mulation [J, |5j are constructed. As is common, our 
derivation starts from considering the evolution equa- 
tions for constraint quantities and looks for sets of data 
for the variables of the main system that guarantee zero 
Dirichlet data for the constraint quantities. The pro- 
cedure is similar to the procedure found in [l2T | but a) 
does not employ reduction to first order, and b) does not 
involve integration of equations in time along the bound- 
ary. Instead, following we rewrite the equations 
in a special form to find well-posed constraint-preserving 
boundary conditions by direct inspection. The approach 
can be generalized to produce boundary conditions of 
the evolving type (see, 0) an d the differential type 
0,mi ■ To further justify the proposed conditions, we de- 
rive an energy estimate for the nonlinear BSS N sy stem 
with boundaries extending the results of 0, Il7| . and 
demonstrate that the nonlinear estimate has the same 
boundary terms as in the linearized case. 

In Section 2 we recall the derivation of the BSSN 
formulation and use the opportunity to discuss choices 
of lapse and shift most commonly found in numerical 
relativity. Section 3 describes the linearization of the 
BSSN equations. In Section 4, the constraint preserving 
boundary conditions are derived, and the generalization 
to evolving and differential boundary conditions is dis- 
cussed. In Sections 5, and 6 the initial-boundary problem 
is defined using the derived conditions. Also, in Section 6, 
a set of boundary conditions for the dynamic part of the 
BSSN system, which is a system first order in time, sec- 
ond order in space, is formulated. Section 7 describes an 
energy estimate for the BSSN system not involving first 
order in space or second order in time reductions. 



II. THE TRACE-FREE DECOMPOSITION OF 
THE ADM SYSTEM 



'Electronic address: alcxander.alcksccnko@csun.edu 
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1 An appro ach not involving first order reduction has been pro- 
posed in llol where boundary conditions were constructed by 
projecting Einstein's equations on timelike boundaries. 



To point out some facts about the nature and proper- 
ties of the BSSN formulation (see 0,0 also, some special 
cases are in 0,^3)' l et us briefly recall the derivation in 
the case of vacuum fields where the right hand side of 
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Einstein's equation is zero. 2 

The derivation starts from the Arnowitt-Deser-Misner 
3 + 1 decomposition 0, |2iJ , 

dohij = -2afaj + 2hi(idj)b l , (1) 
dokij = a[R tj + (kfokij - 2k a k\\ 

+ k u d j b l + ki j d i b l - D.D.a, (2) 

R\ + (A*) 2 - = 0, (3) 

Z>% - Afcj - 0. (4) 

Here a denotes the lapse, the 6j are the components of 
the shift vector b, hij are the components of the spatial 
metric h. The components of the 4-dimensional metric g 
are given by 

goo = -a 2 + bibjh 13 , g oi = b i} g io = h. Lj . 

h lJ denotes the matrix inverse to hij, indices are raised 
and traces taken with respect to the spatial metric; 
<9o := {dt — b s d s ) is the convective derivative, Di is the co- 
variant derivative operator associated to the spatial met- 
ric; the extrinsic curvature kij is defined by equation (I): 
we assume that global Cartesian coordinates t = xq, x\, 
X2, X3 are specified; Rij are the components of the spatial 
Ricci tensor 

R*j — ^f lPg (d p djh iq + didph q j — d p d q hij — didjhp q ) 

~r i<> ti ips 1 - qjr x pqs L ijr )• 

where Tijk are the spatial Christoffel symbols defined by 
Tijfc = (dihkj + djh ik — dkhji)/2. 

The operator in (J2J) contains second order spatial 
derivatives of unknown fields and is very difficult to ana- 
lyze. As a result, it is difficult to judge about properties 
of equation (J2J), and properties of fey in general. How- 
ever, for the trace of the extrinsic curvature k = k\, the 
situation is different. Taking the trace of and using 
0, © we find 

dokj = ak u k H ~ D l Dia. (5) 

The remarkable simplicity of the latter equation suggests 
to separate the evolution of the trace of the extrinsic 
curvature from the system. Specifically, we introduce the 
trace of the extrinsic curvature k = k\ and the trace-free 
part of the extrinsic curvature Aij = kij — (1/Z)hijk as 
new variables. Then (JSJ yields 

d k = \ak 2 + aA lm Ai m - D l Dia. (6) 

Unless the lapse function a is chosen with care, equation 
© is expected to be unstable. For example, for a spa- 
tially independent lapse and zero shift vector, equation 



© yields an estimate dtk > ^ak 2 which implies that 

k > [(1/3) J *a(r)dr + l/fc(0)] _1 , or that the solution k 
is unbounded in a finite time, which is a well-known ex- 
ample of a coordinate singularity. The problem can be 
solved, for exampl e, by imposing maximal slicing in the 
BSSN formulation jig 

D l D ia = ak lm k lm . 

With this condition equation © reduces to dok = 0. 

Alternatively, it is often proposed to use harmonic slic- 
ing 3 which corresponds to setting 

d a = -a 2 k. (7) 

The equation on A is obtained from |3J), J5J, and 

as 

1 2 

d Q Aij = aRtj + -akA i: j - 2aA a A l J + -ak 2 h i: j 
3 y 

- ^aA lm A lm h Z] + A a djb l + Ajidib 1 

-DiDja+^hijD l D ia . (8) 

To proceed with the derivation we need a splitting for 
the spatial metric h compatible to the splitting of k^ 
into k and Aij. In the BSSN formulation, the desired 
splitting is achieved by introducing the conformal fac- 
tor ip = (1/12) ln(det(ft.y)) and the conformal metric 
hij = e~ 4v hij, h 1 ^ — e 4v h l: >. Using Leibnitz formula 
for differentiating the determinant of a matrix 

ddet(hij) = det(h tj )h lm dh ml (9) 

one finds that the derivative of the conformal metric is 
trace- free: 

d~hij = e-^[dhij - ~hijh lm dh lm ]. (10) 

By applying operator do on the definition of tp and using 
0, Q we get the second equation of our system 

d y = -\ak + \d l b l . (11) 
6 6 

Now using HI L)| l and we obtain the third equation 

2 - 

dohij = -2aAij + 2h l{i d j} b l - -h i:j dib l , (12) 

where i„ = e'^A.^, b t = eT^b, are the conformal 
analogs of the variables A and b. Beginning from the 
last equation, indices are lowered and raised with the 



3 Harmonic slicing is a particular case of Bona-Masso family of 
2 A reader not interested in the BSSN derivation may proceed to fe-driving slicing conditions (dt — b l D{)a = —a 2 f(a)k, /(a) > 



Section 3 where the linearized system is listed. l2ll |22| 
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conformal metric hij and its inverse h 1 ^ = e 4v /i u (in this 
case 6 s = b s , and it is easy to redefine d — d t — b s d s ). 

The remaining two equations can be obtained from © 
which can be rewritten in terms of A as 



d Aij = ae~ 4ip R i3 + a(kA i:j - 2A a A l J + -k 2 h 



9" 

1 7/™ 7 7 x 7^7/ 7 ^7/ 2 7 

3' 

3 



e~ Av D l D J a+ e~ 4ip -h l3 D l D ia . 



(13) 



The Ricci tensor in terms of the conformal metric reads 
HE: 



i j ^hP^ {dpdjhiq -\- didphgj dpOqh^j didjhpq^j 



- 2DiDjip - 2h i3 h pq D p D q ip 
_i_ h pq h rs (f f • — f f ■■ ) 

-til it ips-L qjr x pqs^ ijr J 

+ idiipdjip — 4:hijh pq d p ipd q ip. 



(14) 



Here f ljk = (dih kj + d 3 h lk - d k h Ti )/2\ DiVj = diVj - 
h pq Tij p Vq is the covariant derivative associated with the 
conformal metric. The first line in (|14|) can be rewritten 

Rij = ~\h pq d p d q h 3 +d ( Ji pq T ]pq]j) +f pq{l d 3) h pq 

- 2D i D j cp - 2hifh pq b p b q ip + ... (15) 

This suggests to introduce a new variable 

tj = ~h pq T pqj = h pq d p h qj . (16) 

Substituting (|15() in (|13fl one gets the fourth evolution 
equation 

d Aij = -^-ae-^h^dpdgh^ + ae-^d^fj) 



-2ae- 4v DiDjip - 2ae- 4ip h ij h pq D p D q ip 
1 
3 



e-^DiDja + ~ e-^hijrtDta + Wtj, (17) 



where 



Wij = ae-^f pq{z d 3) h pq 

+ ac~ 4<p h pq h rs (r f — f f • 1 

-r it. it \ L i P s L qjr x pqs L zjr ) 

+ 4ae' 4v d l ipd ] ip - 4ae- 4{p h l:j h pq d p tpd q ip 
+ a(kAij — 2AnAj + — k — —A Ai m hij') 

+ A u dfb l + A jl d t ~b l -^A l] d l b l . 

The evolution equation on T is obtained by differenti- 
ating its definition and using the momentum constraint. 
Namely, we apply operator do on (|To|) to get 



5 f z = -2ad l A u -2h pq (d p a)A qi +2aA pq (d p h ql ) 



■ p u,jj-i ql 

+ TA~b l + Ididih 1 + ~h H d s d s b l 



(18) 



Then we notice that h pq D p A iq = d p A pi - T s Af + 
6(d s ip)Af, and thus Q reduces to 

d p A pi - ^dik - f s If + Q{d 3V ,)A\ = 0. 

Solving this equation for d l A\ and substituting the result 
in (|18fl we derive the fifth equation of the BSSN system 



d Ti = --adik + St, 



(19) 



where 



S t = -2aY s Al + l2a(d s ip)At-2h pq {d p a)A q:j 

+ 2aA pq (d p h qj ) + Tidfb 1 + -didib 1 + h h d s d s b l . 

Equations ©, {TU, {El, GZ), an d OH constitute the 
core of the BSSN formulation. These equations are usu- 
ally supplemented by one or more equations describing 
the choice of the gauge functions. Thus, in most cases 
the lapse and the shift are not known but dynamically 
depend on the metric and other quantities. In this work, 
we will assume the harmonic lapse condition (jTJ . Further 
we consider either a prescribed shift b t or a shift that fol- 
lows from the gamma-freezing condition d t Ti = (see, 
for example, (lap. 



III. LINEARIZATION AROUND MINKOWSKI 
SPACE 



Minkowski spacetime in Cartesian coordinates is repre- 
sented by the trivial solution to ADM system: hij = 
kij = 0, a = 1, bi = 0. Consider perturbations of 
ADM variables around the Minkowski spacetime: hij — 
5i 3 + Jij, k i3 = , a = 1 + a, bi = fa, with the 7^, Kij, 
a, and fa supposed to be small. Substituting these ex- 
pressions into the definitions of the BSSN variables and 
neglecting terms of second and higher order in perturba- 
tions we get 



det(hij) = 1 + 7;, <P = j^lL 
0-^ = 1-^1 e 4 ^l + p ; , 



1 



■7ij, 



k — k — \ Aij — A*ij — —SijK, 



r, =: fi = d l %. (20) 

Substituting these quantities in equations (|llfl , l|12H , © 1 
(I17II . i|19|) . and Q and ignoring the terms which are at 
least second order in ip, jij, k, Aij, and Ti we derive the 



4 



linearized BSSN system 

dtV = -\*+\d l Ph (21) 
6 6 

d t a = -«, (22) 

d t n = -d l dia, (23) 

dt% = -2A ij + 2d {i j) -~6 ij d l l , (24) 
d t Aij = -^d l daij + 9 (i r i} - 29i3,-y> 

- 25 ij d l di V - difya + ^5 tJ d l d ia , (25) 

d t T l = -piK+^didfpp + dPdpfa. (26) 

Notice that the linearized harmonic lapse condition is in- 
cluded in this system in the form of equation Q22J1 . In 
fact, this condition will be used in the hyperbolic reduc- 
tion which does not seem possible in general. 

Linearization of Hamiltonian and the momentum con- 
straint equations yields correspondingly 

- 8d l di<p = 0, 9'r z - 8d l diip = 0, (27) 
d l A a - ^d lK = 0. (28) 

Hamiltonian constraint appears in two versions since it 
can be written both in terms of 7 and F. Also, intro- 
ducing the new variable T entails an artificial constraint 

r t = d l j a . (29) 

The linearized problem then consists of determining 
ip, a, k, 7, A, r from equations (|21|I - H26|I provided ini- 
tial data and admissible boundary data. The constraint 
equations l|27|) - (|28|l may or may not be imposed during 
the evolution. The initial data (p(0), k(0), 7(0), A(0), 
r(0) is determined from 7(0) and k(0) using i|20[l . It can 
be checked that if 7(0) and k(0) satisfy the linearized 
Hamiltonian and momentum constraints in the ADM sys- 
tem, 4 then ip(0), k(0), 7(0), A(0), T(0) satisfy the con- 
straint equations (|27|I - H29() . 



IV. CONSTRAINT PRESERVING BOUNDARY 
CONDITIONS 

The BSSN system is a constrained evolution system in 
the sense that it has the dynamic part (|21|l -126 () and the 
constraint part (|27|) - (|29[) . It was assumed for a long time, 
that for the right boundary data, a solution to l|21[) - l|26[) 
will satisfy the constraints automatically once it satisfies 



4 The linearized Hamiltonian and momentum constraints in the 
ADM system are d l d 1 'yii — d l di'y\ = and <9'k« — din\ = 0. 



them initially, however, examples of such data were con- 
structed only recently for a first order reduction. Still 
it remains a question, which we are trying to address in 
this paper, whether a set of constraint preserving bound- 
ary conditions can be proposed for the original BSSN 
variables. 

Let us notice that, for a solution of (|21|I - H26|I . con- 
straints H27fl and (|29|) are consequences of (|28|) . so we 
can focus just on the last one. Indeed, in view of (12 1 1) . 
(|25|l . and (|26l) . the following equations hold for the time 
derivative of (j2T|) 

- 8d l d w ) = d\d l A a - ^diK), 
d t (d l r l -8d l dw) = o. 

These equations state that both parts of (|27|l are satisfied 
as long as they are satisfied initially, and l|28|l is true. 
Similarly, if (|28|l is satisfied, then the time derivative of 
(|2*9"|l is zero in view of J2HJ). Thus (|2*9")l remains zero 
provided it is zero initially. 

We will now construct boundary conditions for system 
(12111 l|26|) that preserve l|28() . We introduce a new vari- 
able 

Mi = d l A a - -diK. (30) 

Equation Ij28(l is satisfied iff Mi = 0, in other words, the 
condition in question must guarantee Mj = 0. 

By differentiating l|30|) twice in time and substituting 
time derivatives of equations (|2~Sf) and ij^ for dfA, <9 t 2 K, 
we derive (terms in tp, a, k, (3, T cancel in view of l|21|) - 

d 2 t M t = d l diMi. (31) 

Initial values M(0) can be determined from A(0), k(0) 
using l|30|l and must be zero for physical initial data. The 
initial values for dtMi can be calculated by differentiating 
(|3l)|l in time and substituting and (|25|) for dtK, and 
dt.Aij , 

^ = 4^^ + ^ + ^-^^. 

It can be verified by substitution, that if 7(0), T(0), ip(0) 
satisfy 123, then d t M{0) l = 0. 

It remains to select the boundary conditions on M that 
imply trivial evolution of Q3ip. However, we notice that 
the boundary data on M is expected not to be given 
freely but determined by the boundary conditions on A 
and k, similar to the way M(0), dtM(0) is determined by 
the main variables A(0), k(0), 7(0), L(0), and ^(0). But 
we do not know how to specify the boundary conditions 
on A and k either! Here is the key: we will select that 
data now by observing its relationship with the bound- 
ary data on M [l3L Hq . For the boundary conditions, 
again, we expect both definition l|30|l and the evolution 
equations (|21|I - H26(I to contribute into the relationship. 
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We introduce scalar products (vi,Ui) = 
Jo ViU 1 dx dy dz and {p^ , <7jj ) = Jn Pi3 a ^ ^ x ^2/ 
for the spaces of vectorfields and matrixfields on ft 
correspondingly. The L2 norms naturally associ- 



ated with these scalar products, \\i> 
\\p\\ 2 = {Pij>Pij)> are denoted by || • || 
energy of system pip. 



\ 2 = (ui,Ui) and 
We introduce the 



e=\\d t M\\ 2 + 110^1 \\ 2 

If we prove that the energy e remains zero at all times, 
then, in view of the trivial initial data, this will imply 
that ||M|| = (e.g, ||). 

Differentiating e in time, using the Green's First Iden- 
tity component-wise to transfer the spatial derivative in 
the second term, and also we obtain 



d,e = 



[ (^-M i )(5 t M i )d<r. 
Jon on 



The energy e is not increasing if 
d 



(^M i){ 8 t M* 



< on dtt. 



(32) 



(33) 



The desired boundary conditions on A and k will follow 
immediately once we rewrite (|33f) in terms of the main 
variables. 

We assume that the boundary 9f2 is a combination of 
arbitrarily oriented planes and consider any of its faces. 
Let vector Hi be the unit vector perpendicular to the 
face, rrii, U complement m to an orthonormal triple (for 
example, m; is any unit vector parallel to the boundary 
and k is the cross product of Uj and m,, k — e{ k njirik — 
(n x m)i). At the flat boundary, the divergence of a 
vector field can be expressed in terms of the directional 
derivatives along vectors n, m, and I, as 

d l v l =n l ^-v l +m l ^-v l + r^-v l . (34) 
an dm dl 

Similarly, the gradient of a scalar field ip reads 

d d d 

d l ip = n.—tf; + rrii—tp + li—ip. (35) 
on am dl 

Next we note that, at any point of the boundary a sym- 
metric trace free matrix is spanned by 

n^nij-j, Tlulj), luvrifj, lilj—mimj, 2ninj — Ijlj — niinij . 

Introducing scalar functions 

Al = 2A ij n^mj), A2 = 2A ij n (i l j} , 
1 



A3 = 2A i n {i m j) , 
1 

~ 6 

we restate A as 



Ah 



AA = 
A 13 (2n l n J 



A 3 



hi j 



A i:i = Al (n^m^) + A2 (n^) + A3 (l^mj-j) 



+ AA (lilj — mimj) + Ah (2ninj 



lilj 



(36) 



liTUj). 

(37) 



Substituting J37|) into J3U|) an d using (gU) and to re- 
place partial derivatives with the directional derivatives, 
we get 



Mi 



l 2dm 



Al 



ld_ 
2dl' 



A2 



2^ Ah- 
on 



2d_ 

3 dn 



A d 1 d ,„ d .. d , c 2fl , 

-— Al + --A3 - — — AA - — Ah - - — « 
2 on 2 ol dm am 3 dm 

\i—A2 + i ±A3 + ^AA - ^Ah - 2 -^n]l- 
[ 2 0n A2+ 2dm Ad+ 0l A4 c^ 5 3 dl^ 

(38) 



The last equation implies that 



(— MA{d t M l ) 

9 A 9 . 1 d ,„ n d <r 2d, 

on 2 om 2 Ol On 3 On 

, r l 9 (1 1 d ,„ „ d .„ 2d, 
x d t {- — Al + --A2 + 2— Ah - -— k] 
2 Om 2 Ol On 3 On 

JLrlA^i + -—A3 _ - —Ah - - — k] 

dn 2 dn 2 dl dm dm 3 dm 

x d [-—Al + -—A3 - — AA ~ —Ah - -—k] 
2 dn 2 dl dm dm 3 dm 

d A d 13 io d AA d Ar 2d, 

+ Tn^Tn A2+ 2^ M+ m M ~dl M -3di K] 



A d 



1 d 



2 d 



x d t \~—A2 + —— — A3 + —AA - — A5 ~ ~—k 



l 2dn 



2 dm 



dl 



dl 



3dl 



(39) 



Either of the two sets of boundary conditions imply 
{(d/d^M^dtMA = on dfl: 



d d 
Al = 0, A2 = 0, —A3 = 0, — AA = 0, 
dn dn 



d_ 

dn 



Ah = 0, 



d_ 

dn 



k = 0, 



dn 



Al = 0, 







A2 = 0, A3 = 0, AA = 0, 



dn 

Ah = 0, k = 



(40) 



(41) 



In particular, l|40|) eliminates the second multiplier in the 
first term of (|39(l and the first multipliers in the second 
and third terms (by commuting partial derivatives and 
using l|46|)'). Condition (|41|) is verified in a similar way. 

More examples of constraint-preserving boundary con- 
ditions can be proposed by inspection of Ij39(l . For exam- 
ple, the condition M^qq = 0, according to JSHJl, is equiv- 
alent to the set of differential boundary conditions that 
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can be implemented numerically P.Hl|: 



1 d 



.41 



2dn 
1 d 



l_d_ 

2dm 
IS 
+ 2dl 



Al 



A3 



ld_ 
2~dl 
d_ 

dm 



A2 



AA 



■2±AS 
on 



_d_ 

dm 



A5 



1 

A2 + -—A3 

2 on 2 dm 



dl' 



-AA 



di 



.45 



2d_ 

3 dn 
2_d_ 
3 dm 
_ 2d_ 
~ 3dl' 



K = 0, 
K = 0, 



0. 

(42) 



Namely, one could prescribe Dirichlet data on ^43, AA, 
and k. Then, l|42[l gives mixed conditions on Al, A2, 
Ab. The problem with this condition is that it is not 
obvious if it leads to a well-posed evolution of l|4t)|) (the 
next example, however, contains an idea on how the well- 
posedness can be established). 

Furthermore, one could have considered a combination 
of Neumann and Dirichlet conditions 



d_ 

dn 



Mirt = 0, Mil 1 = 0, AW = 0. 



(43) 



Applying d/dm to the second, d/dl to the third equa- 
tion of (|42|l and subtracting the results from the normal 
derivative of the first one, we derive (using Ij46|l . I|48|l 
to eliminate {d 2 /dn 2 ) derivatives) an evolution equation 
defined on the boundary 



d 2 



v dl 2 dm 2 
d 2 d 2 
~W dm 2 



d 2 2 d 2 4 d 2 

)A5 ~3W K+ 3 { dl 2 ~ 
d d 



d 2 



dm 2 ' 



)AA 



dl dm 



A3. 



(44) 



Equation (|44ll can be used to find Dirichlet data on Ah 
and k, provided Dirichlet data for A3 and AA is given (see 
[l2|). Once A5 and n are known, the values of (d/dn)Al, 
(d/dn)A2 can be determined from the last two equations 
of (|42(l . The corresponding inhomogeneous algebraic con- 
ditions on Aij then would read 

2{d/dn)A l3 n [l m j) = (d/dn)Al, 
2{d/dn)A ij n [i l j) = {d/dn)A2, 



2A l H 



{l m 3) 



= A3, 



(l/2)A ij (l i l j - mirrij) = AA, 



(l/6)A ij (2n i n j 



lilj 



A5. 



(45) 



The boundary conditions (|45|l are constraint-preservin. 
and represent an analog of conditions introduced in 

El. 



and substituting (|24fl for dtjij in the result, we obtain 
(terms in if, a, k, (3, T cancel in view of (|21l) ~ (|23|l . I|2t)|l ) 



d 2 A tl = tfdtAi 



(46) 



We assume that the initial values A(0) and dtA(0) are 
determined from l|2U|) and l|25|l correspondingly, and that 
either of the conditions (|40|) or l|41|) is given at the domain 
boundary. 

We introduce scalar products (p,, v) — J n [lv da; dy dz 
and (uijk,Vijk) = J^u^kV 1 ^ dxdydz for the spaces of 
scalar fields and triple indexed fields on f). The Li 
norms naturally associated with these scalar products 
are ||^|| 2 = (/x,^) and ||u|| 2 = {u ljk ,u ljk ). The energy 
of system H46fl is defined as 



<* = \{\\9tA\ 



\diAf). 



Similar to Section 4, by differentiating ei in time, inte- 
grating terms with spatial derivatives by parts, and using 
Q46p. we obtain 



d t ex = f (^-Ay)(0 t ^)dcx, 
Jdn dn 



(47) 



Since the right side of (|17|l is zero for either (gUJ or (jlT)) . 

we conclude that t\, and therefore, A remains bounded. 

Similarly, by differentiating l|23|l and substituting 1221) 
for dta one derives the equation for n 



d 2 K 



d l dm. 



(48) 



The boundary conditions (|40|l or (|41|l imply trivial Neu- 
mann and Dirichlet data on k correspondingly. 

Finally, assuming the shift perturbation (3 is known, 
the matrix A is computed from (|46|l , and k is determined 
from (|48l) . variables ip, a, 7, T can be determined from 
(|2Tj) . E2l . pi|l. it^l) by integration in time. 



VI. INITIAL BOUNDARY VALUE PROBLEM 
FOR THE LINEARIZED BSSN. PRESCRIBED 
AND GAMMA-FREEZING SHIFT 



V. EVOLUTION OF A AND k. SECOND 
ORDER IN TIME REDUCTION 

We will argue now that the boundary conditions 14U|) 
and (|41|l lead to a well-posed problem for the linearized 
BSSN system. By differentiating equation (|25|l in time 



Let system l|21ll l|26|) be provided with relevant initial 
data, and (5 be given. If the boundary conditions for A 
and k are taken in either form H4U|) or (|41|) (conditions 
(EESjl and g5j) can be formally imposed and treated simi- 
larly but contain derivatives of the unknown fields), then 
boundary conditions on variables 7, ip, a, and T can be 
obtained by integration of J23|, (|2U, 123, as (the 
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projection operator Nf q is defined below) 



(£3) (or glpj, gg)-|5U) is pen &y 



+ (^]-^)^(°) + /W 
+ - iVf )|-)(29 (p /? 9 ) - j^'A), 
H<p + (1- fi)g-<P = m°) + (1 ~ 



9 9 
fia + (1 - /i)7-a = A*a(0) + (1 - /x)— a(0), 
on art 



(49) 

(50) 
(51) 



nTi =nTi(0)+ / [ 



rT, = t%(0) 



4 l^ + i^ +flP ^: 



(52) 



The projection operator Nf? corresponding to (|40|) is 

Nf q = 2n {p m q) n (l m ]) + 2n {v l q) n (q l j) 

and the one corresponding to (|41|l is 

N pq = 2l i - p m q) l (l m j) + {\/2){l p l q - m p m q ){klj - m, mj ) 
+ (l/6)(2n p n« - Pl q - m p m q )(2n,n j - kl 3 - m.mj), 

/j, = corresponds to (gUJ) and \i = 1 to JSJ. In IpS^. 

Ti stands for vectors Zj, mi. Unless the Neumann data is 
given on n, the equation (|52l) couples T and k through an 
integral equation. If the Neumann data is specified for 
k, for example in (|40l) . the last equation can be replaced 
with 



dn 



4 9 9 
3 9r 9n 



Theorem 1. Let A and k are smooth solutions of 

and corresponding to boundary data {4-Ujl ( or \ 

and the initial data A(0), k(0), dtA(0), 9tK(0) (the last 
two are determined from equations \25\) and ]23\) and the 
initial data 7(0), r(0), a(0), and <p(0) as 



) 



d t Aij{Q) = -\d l da %3 {0) + %r^)(0) - 29;9,^(0) 

- 2S lJ d l d l i P (0) - did ja {Q) + IsytfdtaQ), (53) 
d tK (0) = -d l dia(0). ) (54) 

Then a solution to \21\) - y2b)) satisfying boundary data 



1 1 



^ = ^(0) + / a = a(0) 
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K, 



r, = r<(o) + j + \w p Pp + (55) 

Moreover, if A and k satisfy constraint equation 1286 . 
then 7, T, and if defined by \55)) . satisfy constraints \2 r )\j 
and \2ty as long as they satisfy them at the initial time. 

Proof. Equations J2TJ, US}, and are verified 

by substitution. Replacing 7, T, <p, and a in (|25H and 
(|23|l by their expressions from (|55H and using (I53|) and 
|[5Dl we obtain 

9^ = 9^(0) + / d l diAij, 9 tK = 9 tK (0)+ / 9 ; 9^ 
Jo Jo 

which is a consequence of (O and |@HJl. Substituting 
(|55|l into l|49|l - H52|) we verify the boundary conditions. 

Now consider constraints H27|l . Replacing 7, ip, and T 
with their expressions from l|55|l we obtain 

9*9^- - 8d l dicp = 9'9 J '%(0) - 89*9^(0) 



2&(d l A jl - -djk), 
d l Ti - 8d l d l( p = 9^(0) - 89*9^(0), (56) 



from which it follows that l|27|l is met as long as it is 
satisfied initially and l|28l) is true. Constraint (|29|1 follows 
similarly. □ 

The situation is similar when (3 is to be determined 
from the gamma-freezing condition 9tl\ — which yields 
an elliptic equation for f3 that can be solved at each time 
step 

^d i d p p +dPd p (3 i -^ K = Q. (57) 
The boundary conditions on fi can be taken, for example, 



ftm ! = 0, = 0, — ftn l =0, (58) 



or, 



— /W = 0, ^-^=0, /W = 0. (59) 
on on 

It is beneficial for the computation purposes to replace 
the elliptic equation l|57|) with the hyperbolic equation 




1 
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dtP = ^d p p p + d* d p Pi - - — k, (60) 
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which corresponds to a dynamic gamma-freezing condi- 
tion d t Ti = df(3i. In this case, in addition to H52[) and 
B58J1 , one can consider either of two sets of the radiative 
boundary conditions 

/W = 0, &f = 0, (0 t ft + ^-Pi)n l = 0, 

an 

{d t Pi + fafiiW = 0. (d t j3 i + —(3 i )l i = 0, ^n* = 0. 

After the boundary conditions for ft, A, and k are cho- 
sen, one can set r, — T(0)i = on the boundary, then 
compatible conditions for the rest of the variables follow 
from integration of the system as in the examples above. 
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VII. APPENDIX: ENERGY ESTIMATES FOR 
THE BSSN SYSTEM WITH BOUNDARIES 



for further simplicity zero shift perturbation (Pi = 0), 
the equations J2U, 1|23|) - (|26[) . are restated as 



1 



u 

d t n = —6d l ditp, 

dtjij = -2Aij, 
1 



(62) 

(63) 
(64) 



dtAij = -~d l dt% +d( i Tj ) -88 i djtp, (65) 



(66) 



Taking scalar product of Q63[) with k, integrating the re- 
sult by parts, and using l|62|) to replace k with dt^p, we 
obtain 

ia t [|| K || 2 + 36||^|| 2 ] = -6^(^^. (67) 

Also, we conclude from (EH that 36||<9 t </?|| 2 = ||k|| 2 . 

From (|f]2T) and lj6U)l we observe that d t (Ti — 8ditp) = 0, 
which implies that 

ftllr* - 8d^|| 2 = 0, (d t (Ti - 8di<p),d l %) = 0. (68) 

Next we rewrite the right side of (|65|l in divergence 
form, 



We use differentiation in time to propose boundary 
conditions (|40() , 141|) (and the associated conditions l|49|l - 
(152(1 1 in Sections 4-6. Our derivation, however, strongly 
relies on the linearization assumption and does not ex- 
tend to the nonlinear case. 5 In other words, one could 
have proposed boundary conditions for the nonlinear case 
on the base of (|40|l (or (|41|l ). but one will not be able to 
prove well-posedness of the new conditions by repeating 
the argument of Sections 4-6. In this appendix we try 
to fix this flaw by establishing a well posed energy es- 
timate without reduction to second order in time. We 
use approach proposed by Gundlach and Martin-Garcia 
[l4l Il7j and modify the proof for the case of a bounded 
domain. 

Following 0], we use the densitized lapse 



a = e 6 ^( 



(61) 



which results in significant simplifications in equations. 
The result, however, is expected to extend to harmonic 
slicing and /c-driving slicing as well. 

We illustrate the idea in the linear case, first. Under 
the densitized lapse assumption, equation (|22|l is replaced 
with the linearized densitized lapse condition, the latter 
in the special case Q = 1, reduces to a = 6ip. Assuming 



5 The difficulty appears to be the extra derivatives of the inverse 
metric resulting from differentiation of 1171 that contaminate 
principal part and break the similarity with the linear case. 



d t Aij = d l [--dijij + Si^Tj) - Sdj)?)}, 

then, take scalar product with j4y, and integrate by 
parts: 

\d t \\A\\ 2 + {-\da %3 + 5 l[i {T j) - Sd^ldiAa) 

Replacing A^j with l|64|) in the second term and using 
dt[Siu(Tj\ — 8dj)<p)] = 0, we replace the last identity with 

ld t [\\A\\ 2 + \\^dn l3 ~S l{l (T 3) -8d 3) ^\\ 2 ] 

-l4:% + n (i( r J)- 8d ^)] Aij - (69) 

dn 1 an 

From l|67l) - (|69[) one observes that the energy 
e=|| K || 2 + 36||a^|| 2 + ||r ( -8a^|| 2 + ||A|| 2 

+ \\\da 3l -5 l(l (T 3) -8d 3) ^)\\ 2 , (70) 

has its growth determined by the boundary terms: 

+2 f n {l (T 3) -8d 3) cp)A^. (71) 
Jan 
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Expression (|71(l can be used to propose examples of 
new energy (or constraint) preserving boundary condi- 
tions for the linearized BSSN formulation. However, here 
we will use (|71|l to discuss the meaning of conditions (|4(J|) 
and 1411) proposed in Section 6. Under assumption of 1|40|) 
f (|41|l can be treated similar) and the associated condi- 
tions (|4*9"|l . (|5*Uj) . expression (J7TJ reduces to 

-6— 75(0)A5 + 4(nT,(0) - 8— p(0))A5], (72) 
an an 

where 7 1-5 are the coefficients of decomposition 



metric h), 



(d Q k)ke 4ip = 6 / {d p ip)(d q (ak))h pq - 6 / (fi p d p <p)ak 
6(a p ^)a/c9 9 /i P9 + Ffc]. 



Substituting (|74|1 for ak in the second term and re- 
grouping we obtain our first energy identity: 



ldo[\\k\\l + 36\\d W \\ 2 } = -6 f (n p d p <p)ak + f H 
z Jon Jn 



(78) 



7y = 7I (n (i m i )) + 72 (n(iZj)) + 73 (l^m^) 

+ 74 (Z^Zj — mivrij) + 75 {2niTij — Ijlj — mirrij). 

The energy (|70fl is conserved if the initial data is cho- 
sen as to satisfy (d/dn){p(0) = 0, 712(0) = 0, and 
(<9/<9n)7345(0) = at the boundary (notice that 1^(0) 
can not be given freely but is expected to be subject to 
constraint (|29fl). 

Similarly, conditions (|41|l are constraint and energy 
preserving for $2to~^ if <p = 0, (<9/<9n)7l2(0) = 0, 
and 7345(0) = on dil. 

We are now ready to establish the nonlinear analog 
of l|7Tfl. Substituting (JBTJ into the BSSN equations, dis- 
tributing and expanding covariant derivatives, we rewrite 
the system, 



e 4v d k = -6ad p d q ip + F, 

d a ip = -\ak + -dib 1 , 
6 6 



(73) 
(74) 



dohij = -2aAij + 2h l{l d ]) b l - ^hijdib 1 , (75) 



e^doAij = -ad p [^d p h l3 - h p(l (f f) - 88^)] + G 



doTi = --adik + Si. 



13 > 
(76) 

(77) 



where 



F = 6ah pq f l pq di<p - i8a{d p y){d p <p) - U G ^{d p V )(d p Q) 

-e 6ip (D p D p Q) + \ e 4v ak 2 + e 4v aA pq A pq , 

Gij = -a(d p h p(i )(f j0 - &d j)V ) + Saf p 3 d p ip 
-12a(%>)(6^) + 4ah l3 (d p <p)(d p <p) - 8 {d^d^Q) 

-e^ipibiQ) + I e 6lp h l] (d p ^)(d p Q) 
+ l e ^h i:j (D p D p Q)+ e^Wij. 

Taking scalar product of l|73|l with fc and integrating by 
parts in the right side, we get (n, is the outward nor- 
mal vector to the boundary in the sense of the conformal 



where = J n fc 2 e 4 ^, and 

H = 6(d p <p)[d q d s b s + 6{d q b s ){d sV )]h p i 
+36(d p ip)(d q <p)(aA pq - dH q ^ + \h pq d s l s ) 

+ k 2 e iv {-^ak + ^d s b s ) + &{d p tp)akd q h pq + Fk. 
Also, from (|74() it follows that 

^ |M| 2 = Jj-lak + ^dtb 1 }?. (79) 
Next, we notice from and JZZJl that 



doiTi - 8d iV ) = 8a(di<p)k + I e^(d z Q)k - ^d t d s b s 



-8(dib s )(d s ip) + Si 



Therefore, 



hoWT.-sdMl 2 = I J, (80) 

1 Jn 



J = [8a{d p <p)k + I e 6v {d p Q)k - ^d p d s b s 
-8(d p b s )(d s <p) + S p ](f q - 8d q <p)h pq 
+1^p ~ $d p ip)(f q - 8d qV )[2aA pq - 0<*&«> + h pq d s b s }. 

Finally, taking scalar product of (|76|l with A, integrating 
by parts in the right side, using (|75|l to replace A with 
doh, and re-grouping, we derive our last energy identity: 



[UnJ'dphj) - fi (l {f 3) - 8d j) <p)]aA ij + [ K, 
an z Jn 

(81) 
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where ||i4||2 = /„ A io A ij e 4 ^, and 

K = \\dphij - h p (i(Tj) ~ 8dj)<P)] 
x{[d q (h s ( m d n )b s - -h mn d s b s ) + {d q b s )(d s h mn ) 

-(-2aA q ( m + h lq d( m b l + h s{m d\ q \b s - -h q{m d\ s \b s ){T n) 

4 R 4 
-80 n) <p) - h q{m {8a(d n) ip)k + - e° v (d n) Q)k - -d n) d s b s 

-8(d n fb s )(d s ip) + S n) )]h™~h im h jn 
+aA mn d q Ch pq h m hn} + f -•-;.»''• 



Notice that boundary terms in 1|78|) and i|81|) differ from 
that of linear case in spatial metric only. Notice also, that 
right sides of JJHJ- JEU are combinations of ip, k, h, A, T, 
dip, dh, but not derivatives of k, A, and F, which implies 
that (|78l) - (|ST|l is a closed estimate and may be proposed 
for proving local well-posedness of initial-boundary value 
problem for the BSSN system. In this derivation, we as- 
sumed that covariant components of the conformal shift 
vector, b s , are prescribed, otherwise, terms didjb 8 has to 
be expanded to ensure that they do not contribute to the 
principal part of the equations. 
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